Heat generation and transport due to time-dependent forces 
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We study heat generation and transport properties for solids in the presence of arbitrary time- 
dependent force. Using nonequilibrium Green's function (NEGF) approach we present an exact 
analytical expression of heat current for the linear system. We found that the current can be 
expressed in terms of the displacement of the atoms in the center and the self energy of the heat 
bath. We carry out the calculation for periodic driven force and study the dependence of steady 
state current with frequency and system size for one and two-dimensional systems. We obtain an 
explicit solution of current for one-dimensional linear chain connected with Rubin bath. We found 
that the heat current is related to the density of states of the system and is independent of the bath 
temperature in ballistic transport. The baths can absorb energy only when the external frequency 
lies within the phonon band frequency. We also discuss the effect due to nonlinear interactions in 
the center. 

PACS numbers: 05.40.-a, 44.10.+i, 91.45.Rg, 05.70.Ln 



I. INTRODUCTION 

In recent years, the understanding of heat transport 
in mesoscopic systems has drawn a lot of attention be- 
cause of their interesting physical properties and vast 
applications ranging from nanosize electronic devices to 
thermal transistors. In addition there has been a great 
deal of interest to study how to manipulate and control 
heat. Different theoretical models have been proposed 
to control thermal transport [ll-0|- Several experimental 
works have also been carried out [3, [5| . To understand 
the generic features of these systems numerous studies 
have been done using nonequilibrium Green's functions 
(NEGF) d-Ill, generalized Langevin equation [13, Ell and 
quantum master equation 12, 13] approach. 

The energy transport in general can be achieved as a 
response to the temperature gradient or chemical poten- 
tial gradients and in the linear response regime is gov- 
erned by Fourier's law [13, [3 13 for diffusive systems. 
It is also expected that time-dependent external force can 
induce directed heat transport between the leads at the 
same temperature or even in the presence of temperature 
gradient llq-|l8|. However, whether all energy driven by 
external force can be transmitted to the reservoir or not 
is a valid question. Recent study on driven quantum 
Langevin model for any arbitrary time-dependent poten- 
tial shows that the energy dissipation flow to thermal 
environment is related to the violation of the fluctuation- 
response relation [l^, HOl- Understanding the general 
features of current is therefore one of the main goals in 
nonequilibrium statistical physics. 

For systems driven arbitrarily far from equilibrium it 
is possible to relate the work done during the nonequi- 
librium process with the free energy difference between 
two equilibrium states through Jarzynski's equality (JE) 
l2ll [22i] which states that 



where W is the work done (here the work done is due to 
external time-dependent force) and AF is the difference 
of free energy between final and initial equilibrium states. 
The average is over the work distribution function P{W) 
and /3 = l/(fcsr). If P{W) is Gaussian then it can be 
shown for classical systems that (W) = AF + f3a^/2 
where cr^ = (VF^) — (W)"^ is the variance. An important 
point to realize in this case is that even if the average 
work (W) is independent of temperature the variance in- 
creases linearly with temperature. Finding explicit forms 
of the nonequilibrium distribution functions and hence- 
forth averages for different systems is of obvious interest 
to verify JE 
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In this paper we investigate the influence of the exter- 
nal time-dependent force on a harmonic system which is 
connected with heat baths and analyze the energy current 
with the applied frequency and system size. We explore 
the effect on current due to two different types of heat 
baths, Rubin ^26] and Ohmic ^2^. We discuss briefly that 
one-dimensional linear chain model can not be used as a 
heat pump. To obtain the expression for current and to 
examine the underlying physical process we use NEGF 
method. Aiming for an exact analytical solution of cur- 
rent we consider special form of time-dependent potential 
which is linear in system's position coordinates. 

The paper is organized as follows. In the next section, 
we introduce our model and derive the expression of en- 
ergy current in time domain which in general is true for 
any form of time-dependent force and also in any dimen- 
sion. In Sec. Ill we choose periodic driven force and 
study steady state properties for one-dimensional (ID) 
linear chain and two-dimensional (2D) square lattice. We 
present an explicit solution for current in one-dimensional 
linear chain which is connected to Rubin baths. In Sec. 
IV we discuss the effect on heat current due to nonlin- 
ear interaction in the center. Finally we conclude with a 
short discussion in Sec. V. 
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II. THE MODEL 

We consider an insulating solid where only the vibra- 
tional degrees of freedom plays important role for heat 
transport. Our model consists of a finite harmonic center 
which we denote by C, coupled to two heat baths (L and 
R) kept at temperatures and Tr. For the heat baths 
we consider the standard model of an infinite collection 
of oscillators. Let the displacement from some equilib- 
rium position for the j-th degree of freedom in the region 
a he u'j, a — L,C, R. The Hamiltonian is given by 

'H = nL + nc + nR + nLc + nRc + v{t), (2) 



where 



a^L,C,R 

a — L,R 
(3) 



where superscript T denotes matrix transpose, u" is a 
column vector consisting of all the displacement vari- 
ables in region a, and is the corresponding conju- 
gate momentum. K"' is the spring constant matrix and 
yaC ^ (yCa^T ^ ^) |g ^^le coupling matrix of the 
leads to the central region. V{t) is the time-dependent 
external potential which depends only on the center atom 
variables. In this case the potential has a particular form 
—9{t — to)f^{t)u'-' and f{t) is the time dependent force 
vector acting only on center atoms. The force can be in 
the form of an applied electromagnetic field. For sim- 
plicity we have set all the atomic masses to 1, but the 
formulas can be used for variable masses with a transfor- 
mation Uj — >■ Xj^rrij. We assume that at i < the sys- 
tem is under a known nonequilibrium steady state p with 
respect to the Hamiltonian "Hq where T-Lq is the Hamil- 
tonian without the time-dependent potential V{t). For 
t > to the time-dependent force drives the system into 
a nonequilibrium state. We are interested in calculating 
the current going from the left lead to the center. 

The energy current flowing out of the left lead is given 

by 



hit) = 



dULjt) 
dt 



= -{[n^{t),nm , (4) 



where the average is with respect to the density operator 
p defined above. The operators are in Heisenberg pic- 
ture. The position and momentum operators obey the 
canonical commutation relation 



u'^{t),ul{t)\^ih5jk5''^, a,f3^L,C,R. (5) 



Equation(4) therefore becomes 
hit) = (ulm'^^ucit)) ^ thTr 



CL 



t'=t 
(6) 



/ t 



^2 



FIG. 1: The complex-time contour in the Keldysh formahsm. 
The path of the contour begins at time to, goes to time t, and 
then goes back to time t — to. ti and T2 are complex-time 
variables along the contour. 



Since [uL{t),uc{t)] = 0, the above equation can also be 
written as 



ihTi 



which after symmetrization finally reduces to 
d 



Jr = ihTr 



dt 



-GLcit',t)V 



CL 



(7) 



(8) 



where G(<,t') = ^ [G<{t,t') + G> (t^t')]. The lesser (G<) 
and greater (G^) Green's functions are defined as 



G 



G 



LC.< 
jk 

LC,> 



(<,i') = -^K(i'K(<)) 



jk 



(t,t') 



(9) 



Equation (jS]) is the primary equation we use to calcu- 
late for current flowing out of the left lead. To compute 
the current we need to determine Glc- Our main task 
would be to eliminate the reference to the lead Green's 
functions in terms of the Green's functions of the cen- 
tral region. We use contour-ordered Green's function, 
defined on a Keldysh contour [tI. [28l - [30| (see Fig. 1) from 
to to t and back. The contour ordered Green's function 
can be mapped onto four different normal time Green's 
functions by G""' {t,t') = G{t + iea,t' + iea'), 

where a = ±(1), and G++ = G* is the time ordered 
Green's function, G^^ = G* is the anti-time ordered 
Green's function, G+- = G<, and G-+ = G> . The 
retarded Green's function is given by G''' ^ G*' — G^, 
and the advanced by G° — G^ — G*. These relations 
also hold for the self energy discussed below. It can be 
shown from the equations of motion that the contour or- 
dered Green's function for this model satisfies the equa- 
tion Gcl[t,t') = J dT"Gcc{T,r")V^'^gL{T",T'), where 
the integral is along the contour. The function is the 
contour ordered Green's function for the semi-infinite free 
left lead in equilibrium at temperature T^. Using Lan- 
greth's theorem 30] in Eq. ([5]) we can get 



Ij^{t) = ihTr 



OO Q 

dt" 



to dt 
+ Gcc{tX)mt"'t') 



G^ccit.n^Lit" ~ t') 



(10) 
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with T,L = V'~'^gLV^'-^ being the self energy due to the 
interaction with the left lead. The important point to 
note is that Gcc does not have time-translational in- 
variance because of the presence of time-dependent force 
whereas the surface Green's function g^obeys this invari- 
ance as it is calculated at equilibrium. Our main task now 
is to calculate the center Green's function. 

Let us first consider the one-point contour-ordered 
Green's function for the center which is defined as Q 



F(T') 



Ffr) Ffr^) FJT-) 



(11) 



where Tc is the contour-ordering operator. For one- 
point Green's function contour ordering is not important. 
mP(t) is the operator in the Heisenberg picture. Trans- 
forming to the interaction picture with respect to the 
Hamiltonian T-lo and taking the interaction Hamiltonian 
as V{t) = —9{t — to) f"^ {t)u'-^ we can write the contour 
ordered Green's function as 



Gf(r) 



'Go 



(12) 



where Gq is the Green's function calculated with the 
Hamiltonian Hq . Now if we expand the exponential func- 
tion, the terms with odd numbers of u*^ will be zero since 
the average is with respect to a quadratic Hamiltonian. 
So the expression will contain terms with even number of 
M'^(r) and odd number of /(r) and finally can be written 
in the matrix form as 

G"^(r) = ^ J dT'Go{T,T')f{T') + higher order terms, 

(13) 

(For notational simplicity we have omitted the super- 
script CC on the two-point Green's function of center). 
In Fig. 2 we draw Feynman diagrams for Gp(r) upto 
third order of force. The contribution from the first di- 
agram is nonzero. However, all the higher order terms 
contain the same type of vacuum diagrams which are 
zero. Vacuum diagram in this case is defined as a dia- 
gram where all variables are integrated and the result is 
independent of space or time. The expression for such a 
diagram in terms of contour variable can be written as 



dTdT'fiT)G^^{T,T')f{T') 

I <ydta'dt'ni)''Gr'Mr\t'). (14) 



The last line is obtained by going to the real time using 
Langreth's rule. Since the driven force / does not depend 
on the branch index f'^{t) = f~{t)— f{t), we can take 
the summation inside and obtain [7| 



aa'G''/ - G* + G* - G< - G> = 0. (15) 

So the above expression is zero. Similarly all the 
higher order diagrams doesn't contribute to the one-point 



i,T j,T' i,T j,T' I,T" m,T"' 



higher order terms 



FIG. 2: The Feynman diagram for one-point Green's function 
of the center in the presence of time-dependent force. 



Green's function. So the exact expression for (uc(t)) is 
now given by 



{uc{t)) 



dr'Go{r,T')f{T'). 



(16) 



From this expression it is also clear that {uc{t)) does not 
depend on the branch index i.e. (u^(<)) = (wp(i)). So 
in real time we obtain 



{uc{t)) = - dt'Gl{t-t')f{t') 



(17) 



where Gq is the retarded Green's function and is defined 



Gl,^{t.t')^--e{t-t'){[u^(t),u'i{t')]). (18) 

It is also related to the response function in the linear 
response theory. In fact, the same result, Eq. (17), can 
also be derived from the standard linear response theory. 
Similarly the two-point Green's function in the interac- 
tion picture is also calculated using the definition and is 
given by 



G,.(r,r') = -A(^,^.f(r)u^(r')e5:™i/'^^"/^(^'>™(^"))^^. 

(19) 

As discussed above we can expand the exponential and 
the terms greater then 0{p) vanishes as they contain 
same type of vacuum diagrams. The exact expression 
can be written as 

G-,fe(T, r') = G'ojfe(T,T') - ^ / dTldT2Gojm{T,Tl) 

Go,fes(r',r2)/„,(Ti)/,(r2). (20) 

In terms of {uc{t)) the center Green's function now be- 
come 



G(r,r') = Go(r,T') - -{uc{T)){uc{T')f . (21) 

From the above equation we can write G — Gq + 6G with 
SG — —j^{uc{t)){uc{t'))'^- Now using the property of 
{uc{t)) we can write 



5G++ = 5G+- ^ 5G-+ = 5G- 



(22) 



which implies that 5G'' = SG" = Q and 5G< = 6G> = 
5G = -{{uc{t)){uc{t'))'^ . So using Eq. (10) the expres- 
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sion for the current reduces to 



I Lit) = ihTr 



dt' 



dr 



to 



tL(t' -t")G-^{t\t) 



+ TTdt' -t")G^{t",t) 



+ Tr 

= im+iiit) 



dt"{uc{t)){uc{t")r^mt"^t') 



(23) 
t'=t 



By writing Ihit) in this form it is clear that the contri- 
bution to the energy current is separated into two parts. 
if^it) is the current due to driven force and /£(i) is due 
to the temperature difference between the heat baths. In 
the long time limit i.e, t oo, If^{t) is the steady-state 
heat flux and is given by the Landauer like formula 0| 



11 = 1- 



dujhujT[uj] [Jl ~ Ir) 



(24) 



where fa = 1/(5'^°''" — 1) is the Bose-Einstein distribu- 
tion function where Pa — ^/{ksTa), T[w] is known as 
the transmission function and is given by the Caroli for- 
mula TH = Tr[G5rtGgr/j] with Ta ='i(Sj; - Sa) and 
Gq[uj] = [Go[w]]\ a = L.R. The separation of energy 
current into two parts is possible because the system is 
linear and the driving force is not correlated with the 
heat baths. 

If we take the two heat baths to be at the same tem- 
perature, i.e, AT = Tl — Tr = 0, then /£ is zero. So 
in the linear case the final expression for current with 
AT = is 



/^(t)=Tr 



'dt"{uc{t)){uc{t")r-l;^l{t",t') 



J t'=t 
(25) 

where i;^(i" - t') = if t" - t' > 0. This is the central 
equation which can be used to calculate the current both 
in transient state as well as in steady state with arbitrary 
form of force. This expression is true for systems with 
finite heat baths and also in higher dimensions. 

In the following we will consider situation for AT = 
and use Eq. (^5]) to calculate the current. We take a 
particular form of force which is oscillatory and carry 
out calculation for ID linear chain and 2D square lattice 
for two types of heat baths (1) Rubin bath and (2) Ohmic 
bath. 



III. PERIODIC DRIVEN FORCE 

We consider the form of force given by f{t) = /oe^*^^* + 
c.c where /q is a column vector with complex ampli- 
tude and is the driven frequency. Then from Eq. (jl7p 
{uc(t)) can be written as 



where GJjfri] is given by 

Glin] = [{n + ij^fi -K^~ - ^R[n]\ (27) 

with 7] — > 0+ and / is the identity matrix. We set 
to —00 for steady state oscillation and finally average 
over a time period II = ^ /J" lL{t) dt where t = 27r/fi is 
the time period of the driving field, we finally get from 
Eq. dlSl), 

II = -nS[n], (28) 
S[n] = Tr(G5[fl]/o/tGSPri[l]]), (29) 

Since ^^^[ri] is always positive the current is fiowing into 
the lead. The average rate of work done is positive and 
consistent with the second law of thermodynamics. For 
this particular case the same result can also be obtained 
using linear response theory. We can write Eq. (28) in 
another form by using the following relation between Gg 
and Gg 

GS[0] - G^oM = -iG'op] {FLM+TRin]) G^p] (30) 



then we can write 



(31) 



which is a consequence of energy conservation and Ic — 
zOTr[(G5H-Gg[C!])/o4]. 

It is important to note that the above expression 
(Eq. (29)) contains Gg, which are independent of tem- 
perature. So in the ballistic case the current is indepen- 
dent of the temperature of the heat bath. However, the 
higher moments of current, for example (/£), in general 
do depend on temperature. 

For the Hamiltonian given in Eq. (2) it is also possible 
to calculate work done by the external time-dependent 
force which is given hyW = - /J" dtj'^{t)uc{t) where the 
dot refers to derivative with respect to time. Using this 
definition one can then calculate (W) and (yV^) and can 
verify JE Following this definition, P{W) is Gaus- 
sian and equivalent statement of JE classically reduces 
to {W) = § {{W^) — (W^)^] . Since the integration is over 
a time period r = 2^/0, the initial and final equilibrium 
states are the same and hence AF = 0. It is also impor- 
tant to realize that if we define W = - /J" /i,(i) dt then 
it does not satisfy Jarzynski equality and P{W') is not 
Gaussian. However the relation between first and second 
moment come out to be the same classically. 

The first and second moment of W (only the driven 
force contribution) can be written down explicitly 

{W) = T^S[^V\ 



(fLm^fRm) 



(32) 



(ucit)) = Gl,[n]foe 



C.C. 



(26) 



When the leads are at the same temperature, we have, 

(Vt/'2) _ {wT = fi^il + 2fLm) (W), (33) 
which classically reduces to (W^) - {W'f = § (W). 
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FIG. 3: Energy current II as a function applied frequency 
for different system size of one-dimensional chain with force 
/iW = i-iy foe-'""' + C.C. (a) iVc=4, (b) Nc=6, (c) Nc^8, 
(d) Afc=10. K=l eV/(uA2). 
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FIG. 4: Energy current II as a function applied frequency 
for different system sizes of one-dimensional linear chain with 
force fj{t) = i-iy foe-'"' + c.c. (a) Nc=l, (b) Nc=3, (c) 
Nc=5, (d) Nc=7. K=l eV/iuK^). 



A. Application to ID chain 

1. Rubin bath 

Here we consider a ID chain with inter-particle spring 
constant K. We divide the full infinite system into three 
parts, the center, the left and the right lead. The leads 
are at the same temperature with the center. We drive 
the center with the force f{t) and evaluate Eq. (^5]) . The 
classical equation of motion for the center atoms is given 
by 



K{uj^i - 2uj + + fj{t), l<j<Nc (34) 

where Nc is the number of particles in the center. The 
leads obey similar equations with fj{t) = 0. The equi- 
librium Green's functions satisfy time-translational in- 
variance and hence Fourier's transform exists. In fre- 
quency space the retarded Green's function for the semi- 
infinite linear chain can be obtained by solving [31 1 
[(il + iiff' — K]G''q = /, where matrix K which is infi- 
nite in both directions is 2K on the diagonal and —K 
on the first off-diagonals. The solution is translationally 
invariant in space index and is given by 



\\j-k\ 



(35) 



with A 



AK"^ and u) = {^l + i-qY 



2K, 

Choosing between plus and minus sign by |A| < 1. The 
surface Green's function in frequency space is given by 
S'£[i7] ~ ~KX. It is clear from the expression of A that 



2K 2K 



it is complex within the range < fl < 2^/ K and is real 
outside this range. Hence Fl is zero outside the phonon 
band. 

Here we consider the force fj{t) = /^e~*^* -|-c.c. where 
= (— l)-'/o which also mimic the structure of a crystal 




FIG. 5: Energy current II versus length of the center for 
different applied frequencies for one-dimensional linear chain. 
Here (a) J7=0.39, (b) 0=0.78, (c) 0=0.98, (d) 0=1.40. The 
frequencies are given in 10^* (Hz) unit. The other parameters 
same as in Fig. 3 . 



having alternate charges at the sites. For this force the 
expression for current is 



(l-(-l)"c cos(Ncq)) 



II = 



2K 



sin q ( 1+cos q) 



.0, 



for < f7 < 2^/K, 
for n > 2V^. 



(36) 

II is of order 1 and q is given by the dispersion relation 
OF- = 2K{1- cosq). 

In Fig. 3 and Fig. 4, we plot energy current as a func- 
tion of applied frequency for different system size. The 
value of force constant is chosen as if = 1 eV/(uA^) and 
/o = 1 nN in all our calculation. In Fig. 4, the current 
is nonzero even at zero frequency because the system as 
a whole is not charge neutral. For Nc = 1 the cur- 
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rent II is proportional to the density of states (DOS). 
More importantly the current is exactly zero when the 
applied frequency matches with the normal mode fre- 
quency of the system and the corresponding wave number 
is given by for even Nc, q = 21:111/ Nq and for odd Nq^ 
q = (2n+l)7r/A^c withn = 0, 1, • Therefore the 
number of resonance peaks and number of zero's depends 
on the eigenfrequencies and hence on the size of the cen- 
ter system. The average current diverges at 51 = 2\/K 
as the DOS of the full system diverges at the maximum 
frequency of the whole system. For fl > 2\[K the system 
does not allow energy to pass through. Similarly one can 
calculate the right lead current 7^ and the expression is 
the same with Eq. ([36]) . Since we apply force on all the 
atoms of the center by symmetry argument we can say 
that the total input current Iq divides into two equal 
parts and goes into the leads i.e. = \Ib\ — |/c|/2. 

In Fig. 5, we give results for energy current as a func- 
tion of total number of particles in the center for differ- 
ent values of external frequency. For finite systems the 
current oscillates with system size and depending on the 
values of $7 it shows periodicity with respect to Nc- The 
maximum amplitude of the average current is fixed and 
is proportional to VLf^/2K. 



2. Ohmic bath 

Here we consider the center system to be connected 
with two Ohmic baths. The difference between Rubin 
and Ohmic bath is that, the self energy in this case is 
approximated as = ijQ where 7 is the friction coef- 
ficient. More precisely the Sl and S/j matrices are given 
by 



i^fQSijn Sin- 



(37) 



Using the form of the Green's function in Eq. (27) and 
after some bit of algebraic simplifications, we obtain the 
following results 



where 



lL^2jn'f^\g[ 



Nc 

5P=^(-l)^+^G^,-[f7] 



(38) 



(39) 



From the above expression and Eq. (27) it is clear that 
energy current depends on the denominator A[r2] = 
|det[D[0]]|2 where£i[f7] = (Q^ I- K'^ -inVL-inT r) is 
Nc X Nc matrix. The matrix elements are given by Dij ~ 
Sij{n'^ -2K - in"f{6,,i + 5i^M)) - K6^,j+i - K 5,.j-i- 
If we denote PAr£;[fi] = det(57^ — K'^) to be the charac- 
teristic polynomial of the matrix K'-^ with Nc particles 
then it can be shown that [32*1 



A[n] = [PncM 



7' ^''Pnc- 



47^ pf 



Nc 



(40) 
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FIG. 6: Energy current II as a function applied frequency 
for different values of friction coefficient 7 of one-dimensional 
linear chain with force fj{t) = (— 1)^ /oe"'"' ~f c.c. (a) 7=0.01, 
(b) 7=0.5, (c) 7=3.0, (d) 7=5.0. K=l eVI{uk^) and Nc =8. 



where Pnc-i[^] is the polynomial of the (A^c — 1) x {Nc — 
1) force constant matrix K'-^ with first row and column 
or last row and column taken out from K'~' and similarly 
PNC-2M is the polynomial of the {Nc - 2) x {Nc - 2) 
matrix by taking out the first and last rows and columns 
from - The resonance and the zero's of current cor- 
responds to the minimum and maximum value of A[Q] 
respectively. It is difficult to obtain explicit solution in 
this case. However the equation become simple for small 
and large value of 7, the friction coefficient. For small 
friction it is clear from Eq. (40) that A[Q,] = P^^[il]- So 
the resonant frequencies depends on Nc eigenfrequencies 
of the force constant matrix K'^ - In the opposite limit 
i.e, for large 7 we obtain A[i}] = Pnc-2[^]- depend- 
ing on the value of 7 the resonance peaks shift from Nc 
to Nc - 2- 



In Fig. 6, we plot the current with applied frequency 
for different values of damping coefficient 7. The value of 
7 is chosen in proper units. The zero values of the current 
is same as in Rubin's case. However there is a gradual 
shift in the resonance peak depending on the parameter 
7. The current doesn't diverge at ri = 2y/K and the 
width of the peaks depends of 7. We check numerically 
the behavior of 1^ with system length and we found that 
the behavior is similar with Rubin baths. In this case 
also we have = \Ir\ = \Ic\/2- 

Similar Ohmic model was also investigated by Marathe 
et. al [33! for Nc — 2 where they conclude that this model 
cannot work either as a heat pump or as a heat engine. 
Our calculation agrees with their results. 

It is also possible to calculate current in the over- 
damped regime by dropping the term (fi -I- irj)'^ in ^[11] 
given in Eq. (27). In this regime for N — I our result 
agrees with the result obtained in Ref. [2^ for magnetic 
field B = 0- 
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3. Comparison between Rubin and Ohmic bath for driving 
force on single site 

As we have seen that if we apply force on all the atoms 
of the center because of the symmetry of the problem 
if we interchange the left and right lead (which we as- 
sume to be the same) the value of the current should 
not change and hence we have the only possible solution 
— \Ir\ = But this is not the case, at least 

for Ohmic bath, if we apply force on a single or multi- 
particles but not on all. 

If we consider the force on the ath particle as fi{t) = 
dia (/q e~*^* -I- c.c) then for the Rubin bath case using 
Eq. (28) and Eq. (29) we get 




0.5 1 1.5 

0(10''' Hz) 



0.5 1 1.5 2 

n(io"' Hz) 



h - -2r!XIm(A)/o"(/o")*|Go,oi 
In = -2nKlm{X)f^{f^y\Go.Nc 



(41) 



Using the solution for G'Q[il] given in Eq. (35) we obtain 



FIG. 7: Energy current II and Ir as a function of applied 
frequency for driven force at different site of one-dimensional 
chain connected to Ohmic bath, (a) and (b) are for a—1 and 
(c) and (d) are for a = 3, Nc =16. K = 1 eV/(uA^). 



n 



2K sin q 



JSif^y, forO<r!<2Vx, 
[0, for > 2y/K. 

(42) 

which says that, because the full system is translationally 
invariant in space, the magnitude of current does not 
depend on which site the force is applied and hence | /l | = 
\Ib\ = l-^cl/2 is the only possible solution. The result is 
similar with Nc = I in Eq. (I36p . 

However, this scenario is not valid for Ohmic bath. In 
this case the full translational symmetry is broken and 
hence applying force on different sites generate different 
magnitudes of current on left and right lead. In Fig. 7, we 
plot the heat current 1^ and Ir for one-dimensional chain 
as a function applied driving frequency at different sites. 
Clearly 7^ and /r are different in magnitudes. Hence by 
applying force on different sites it is possible to control 
current in both the leads for Ohmic case. 



4- Heat pump 

Heat pump by definition transfers heat from cooler 
region to hotter region. One-dimensional linear system 
with force applying on any number of sites fails to work 
as a heat pump. To understand the reasoning we may 
consider Eq. (23) which says that heat current II is a 
sum of two parts. If we assume > Tr then the first 
term in Eq. (23) which gives the steady state heat flux 
due to temperature difference is positive, i.e, current goes 
from left to right lead and the driving term which does 
not depend on temperature, always contribute a negative 
value to both I^ and Jr. Hence Ir is always negative in- 
dependent of whether we apply force on one site or on all 
the sites. So it is not possible to transfer heat from right 
lead to left lead in this case. 



B. Application to 2D square lattice : Rubin Bath 

In this case we consider a square lattice with force con- 
stant K both in x and y direction. We take a small part 
of the full infinite system which is square in shape and 
call it the center and rest is treated as a bath and is kept 
at a constant temperature with the center. The classi- 
cal equation of motion for the x-component of the center 
atoms is given by 



l,k 



l<j,k<Nc; 



(43) 



and similar equation for the y-component. The total 
number of particles in the center is N^. The retarded 
Green's function for the full system is given by [3^ 



1 



^ gi(Ri-R,,)-k 

^"'[^] = N^1^Q2_2K[2- cos(fc^a) - cos{kya)] ' 

(44) 

where I = I2 + {h - 1)N and I' ^ I2 + {I'l - l)N, R; = 
Ziai + hsi-i, R;' = I'lSii + ^332 and N is the total number 
of particle in the full system, k is the wave- vector and it's 
components are given by — ^j^i ky — '^JPj where a 
is the lattice constant, ai and a2 are the primitive lattice 
vectors. In the large N limit, i.e., N — > 00, one can write 



1 

2^ 



dqy cos{n2qy) L{qy), (45) 



where ni = l[ — h, n2 — 1'2 — h and 



27r 7-7r ^ 2_ft' 1 2 — cos(qx) — cos{qy 



(46) 
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1 2 

acio" Hz) 



0.5 1 1.5 2 2.5 3 
£2(10'" Hz) 



FIG. 8: Energy current / as a function applied frequency 
for different system size of square lattice with force Fj(t) = 
(-l)jfoe-^"* + c.c. (a) jVc=4, (b) Nc=Q, (c) iVc=8, (d) 
Nc=lQ. K = 1 eV/(uA^). 



This integral can be done using the contour integration 
technique and can be written as 



L{qy) = 



(47) 



with A 



± 



and uj 



2K (2 



2K 2K 

cos(gy)). The choice between plus and minus sign de- 
pends on |A| < 1. In this case the explicit expression 
for S2"[i^] is not known. However, by knowing G[;/[r2] 
we can compute the self energy of the infinite 2D square 
lattice with a removed square part using the following 
equation 



lll[^] = {^ + i^f-Kc-[Gl[n]\ 



(48) 



We can obtain Tl from the above expression as ri:,[0] = 
-2Im(E£[f2]) = 2Im[G'5[ri]]"\ 

In Fig. 8 we plot the average current going out of the 
center with frequency where the force is in both in x and 
y direction with same magnitude (/o = 1 nN) and is given 
by fj(t) = (— l)jfoe~'^* + c.c. The behavior of average 
current in this case is quite similar to the ID case. The 
oscillation also increases with Nc and the value of current 
goes to minimum when the applied frequency matches 
with the normal mode frequencies of the full system. 

In Fig. 9 we plot the current with system size and it 
is found that the current oscillates with Nc and it also 
shows a periodic pattern depending on the value of fi. 
Tl is roughly proportional to N^, the total number of 
particles in the center. 



IV. NONLINEAR INTERACTION AT THE 
CENTER 

It is possible to study the effect due to nonlinear inter- 
action in the center for this model. In this case we assume 




FIG. 9: Energy current I versus size of the center for different 
applied frequencies for two-dimensional square lattice. Here 
(a) t7=0.10, (b) Q.=IA7, (c) !^=0.98, (d) 17=0.39. The fre- 
quencies are given in units of 10^'' (Hz). The other parameters 
same as in Fig. 3. 



that both the force and the cubic interaction switched on 
at f = — 00. So using contour-ordered Green's function 
and interaction picture the center Green's function can 
be written as 



(49) 



We take cubic potential which is 

Hn{r) = \Y.j j dT"Tjkl{T,T',T")Ujir)UkiT')uiiT") 
jkl 

(50) 

where Tjki{T,T' ,t") = TjkiS{T,T')6{T,T"). If we expand 
the nonlinear interaction the first term gives us our old 
linear result. The first nonzero contribution comes from 
the second term of both the nonlinear potential and the 
force. The expression for the Green's function in first 
order of force is given by 



J2iTrUm{T)ur,{r')fo{r")uo{r") 

jklo 

Tj-fei (ri , T2 , T3 ) (ri ) life (t2 ) (t3 ) ) Go ■ 



(51) 



Since the density operator is quadratic we can use Wick's 
theorem and finally wc get 15 possible terms which gives 
rise to three independent Feynman diagrams. The final 
expression combining all these diagrams is 



-I fOO /"OO 

lL{t) = -^ E ^i*^' / "^"^l duj"h^'e-'^"'G^^^[J'] 

jklm.no 

Gkn[u^']T.l^[uj'] + GmiW +uj"]GUJ\>:i^[J]) (52) 
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From this expression it is clear that in the steady state 
there is no contribution to current to the hnear order 
in / if we consider a cubic inter atomic potential. To 
see the effect due to nonlinearity and also temperature 
dependent heat current it is important to go to higher 
order in force and also of the nonlinear potential. 

V. CONCLUSION 

In summary, we present an exact analytical expression 
of energy current for driven linear system in time do- 
main. The energy current is written in terms of the dis- 
placement of the center atoms and self energy of the heat 
bath. We study the properties of energy current for two 
different types of heat baths with different forms of self 
energy E. We obtain an explicit expression of current for 
one-dimensional linear chain, connected to Rubin baths, 
exploring the translational symmetry of the full system. 



We discuss the similarities and differences between Ru- 
bin and Ohmic bath when the force is applied on all sites 
or on single site. We also relate the time integral of left 
lead current with work and discuss that this particular 
definition does not obey JE. However, we find that the 
relation between first and second moment of work in both 
cases are same, classically. It will be interesting to study 
the general features of current using Eq. (25) with other 
forms of time dependent forces. The effect on current 
and heat pumping due to nonlinear interaction in higher 
order of force are worthy of further explorations. 
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